function [output] = fn_s_hat(R1,lambda,s_min,s_max,phi)

% myfun = @(s) (R1-fn_rho1(s,phi))./(1-1./fn_rho2(s))-(1-lambda)*R1;
% 
% x0_min = 1.00001;
% x0_max = 1000*s_max;
% x0     = [x0_min, x0_max];
% 
% if (myfun(x0_min)>0 && myfun(x0_max)>0)
%     output = s_max;
% elseif (myfun(x0_min)<0 && myfun(x0_max)<0)
%     output = s_min;
% else
%     options      = optimset('Display','off');
%     [sol,~,flag] = fzero(@(s) myfun(s),x0,options);
%     output       = min(sol,s_max);
%     if flag ~= 1; disp('Error solving s_hat'); end
% end

zeta = 1;

C(1)  = phi;
C(2)  = -(R1*(1-(1-lambda)*zeta)-1+phi);
C(3)  = -(1-lambda)*R1*zeta;
sol   = roots(C);
s_sol = max(sol);

output = s_sol*(s_sol<s_max && s_sol>s_min) + s_min*(s_sol<=s_min) + s_max*(s_sol>=s_max);

% zeta = 1;
% 
% C(1)  = phi(2);
% C(2)  = -(R1*(1-(1-lambda)*zeta)-phi(1));
% C(3)  = -(1-lambda)*R1*zeta;
% sol   = roots(C);
% s_sol = max(sol);
% 
% output = s_sol*(s_sol<s_max && s_sol>s_min) + s_min*(s_sol<=s_min) + s_max*(s_sol>=s_max);

end